R Notebook

Data from 1960 to 2010 raining days, there are 4168 observations 13 days of which are rain exceeds 1000 (10cm)

Change of exceeds 1000 is 0.003119002 0.31% change

library(tidyverse)
library(readr)
library(drc)
library(reshape2)

#install.packages("moments")                
library("moments") # for skewness of data

Data

read.csv("./raw/prcp.csv") %>%
  mutate(date = as.Date(date)) %>%
  mutate(year = format(date, "%Y")) %>%
  filter(prcp != 0) %>%
  dplyr::select(date, year, prcp) -> data

data %>%
  filter(year <= 2020) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist

Modling Distribution

Use data from 1960 - 2010 for modeling precipitation

hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 

data.frame(x, y) -> dens 

summary(model)

Model fitted: Exponential decay with lower limit at 0 (2 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
d:(Intercept) 0.1407885  0.0012481  112.80 < 2.2e-16 ***
e:(Intercept) 1.3292414  0.0077504  171.51 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.001397153 (945 degrees of freedom)

plot(x = dens$x, y = fitted(model), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

mod.exd <- drm(y ~ x, fct = EXD.2(), data = dens)

plot(x = dens$x, y = fitted(mod.exd), type = "l", col = "blue", xlim = c(0,300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

summary(mod.exd)

Model fitted: Exponential decay with lower limit at 0 (2 parms)

Parameter estimates:

               Estimate Std. Error t-value   p-value    
d:(Intercept) 0.1167687  0.0021435  54.476 < 2.2e-16 ***
e:(Intercept) 9.3269006  0.2794558  33.375 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error:

 0.00182667 (945 degrees of freedom)

1 - pexp(log(1000)*k, d)

[1] 0.005534439

1/d

[1] 7.102853

data %>%
  filter(year <= 2010) %>%
  pull(prcp) %>%
  mean() %>%
  log()*k

[1] 23.28785

0.31% change versus 0.442%

Labs:

A Single Years 20 year roll

yr = 1960
data %>%
  filter((year >= yr) & (year <= yr + 9)) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 
data.frame(x, y) -> dens 
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)

[1] 0.0002201357

plot(x = dens$x, y = fitted(mod), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

yr = 1970
data %>%
  filter((year >= yr) & (year <= yr + 9)) %>%
  pull(prcp) %>%
  hist(breaks = 1000, plot = FALSE) -> hist
hist$counts %>% sum() -> sum
x <- hist$mids 
y <- hist$counts/sum 
data.frame(x, y) -> dens 
mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
d <- mod$parmMat[1]
e <- mod$parmMat[2]
k <- 1/e/d
1 - pexp(log(1000)*k, d)

[1] 0.002985969

plot(x = dens$x, y = fitted(mod), type = "l", col = "blue", xlim = c(0, 300))
lines(x = dens$x, y = dens$y, type = "p", col = alpha("black", 0.5))

1 - pexp(log(1000)*k, d)

[1] 0.00457532

mean = 76.84094 sd = 144.0488

Explonential Explainability

Risk Trends

Risks Through out Years

mtr <- data.frame(year = YEAR)
list <- c()

for (i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 9)) %>%
      pull(prcp) %>%
      hist(breaks = 400, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$counts/sum 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(1000)*k, d)
    list <- c(list, p)
}

NaNs producedNaNs producedNaNs produced

plot(1:length(list), list)

Use smaller breaks

Find the right break numbers

YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22)

lenlist <- c()
kulist <- c()
for(yr in YEAR){
data %>%
      filter((year >= yr) & (year <= yr + 19)) %>%
      dplyr::select(year, prcp) %>% 
      pull(prcp) %>%
      length() -> n
  lenlist <- c(lenlist, n)
data %>%
      filter((year >= yr) & (year <= yr + 19)) %>%
      dplyr::select(year, prcp) %>% 
      pull(prcp) %>%
      kurtosis() -> ku
  kulist <- c(kulist, ku)
}
Square_root = c()
Sturges = c()
Rice = c()
Donna = c()
for(n in lenlist){
  sq = ceiling(n^0.5)
  sg = ceiling(log(n, 2)) + 1
  rc = ceiling(2*n^(1/3))
  igh = (6*(n - 2)/(n +1)/(n +3))^0.5
  dn = log(abs(1 + kulist[which(lenlist %in% n)]/igh) , 2) + log(n, 2) + 1
Square_root = c(Square_root, sq)
Sturges = c(Sturges, sg)
Rice = c(Rice, rc)
Donna = c(Donna, dn)
}
data.frame( Year = YEAR,
          Length = lenlist,
           Square_root = Square_root, 
           Sturges = Sturges,
           Rice = Rice, 
           Donna = Donna
           ) -> meth

meth

Work Flow

In a 10 years interval compute risks of precipitation higher than 10cm (prcp > 1000) by using different breaks of histogram from usgin 250 to 1000

Choosing bin have many choices actually

YEAR <- c(1960, 1980, 2000)
BREAKS <- c(19, 20, 21, 22, 30, 40, 50, 60, 70, 80, 90)

for(n in BREAKS)
  data %>%
      filter((year >= 1960) & (year <= 1969)) %>%
      pull(prcp) %>%
      hist(breaks = n)

kurtosis(data$prcp)

[1] 31.3595

list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
klist <- c()
for(n in BREAKS){
  for(i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 19)) %>%
      pull(prcp) %>%
      hist(breaks = n, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$density 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(1500)*k, d)
    list <- c(list, p)
    listB <- c(listB, n)
    Year <- c(Year, i)
    dlist <- c(dlist, d)
    elist <- c(elist, e)
    klist <- c(klist, k)
  }}
mtr <- data.frame(year = Year,
                    Breaks = listB, 
                      Probability = list, 
                  parameter_d = dlist,
                  parameter_e = elist,
                  parameter_k = klist)

Visulisation

ggplotly(g)

n too large, allowed maximum for palette Blues is 9
Returning the palette you asked for with that many colors

mtr %>%
  mutate(year = as.factor(year)) %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(y = parameter_d, x = parameter_e, color = year)) + 
    geom_jitter(alpha = 0.8) + 
    scale_color_brewer(palette = "Reds") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        ) + 
    xlab("parameter e: skewness of the curve") + 
    ylab("parameter d: height of the curve") 
    

Parameter d influence the height of pdf, parameter e influence the lenghth of pdf… this suggest in general, precipitation distribution cure are getting flatter and

plot_ly(x= mtr$parameter_e, y= mtr$parameter_d, z=mtr$Breaks, type="scatter3d", mode="markers", color=mtr$year)

mtr %>%
  dplyr::select(year, Breaks, parameter_d, parameter_e) %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  melt(id = c("year", "Breaks"), variable.name = "parameter") %>%
  ggplot(aes(x = year, y = value)) + 
    geom_line(aes(x = year, y = value, linetype = parameter, color = Breaks)) +
    scale_color_brewer(palette = "Purples") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

mtr %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(x = year, y = parameter_k, color = Breaks)) + 
    geom_point() + 
    geom_line() + 
    scale_color_brewer(palette = "Blues") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

calculate different

list <- c()
listB <- c()
Year <- c()
dlist <- c()
elist <- c()
for(n in BREAKS){
  for(i in YEAR){
    data %>%
      filter((year >= i) & (year <= i + 9)) %>%
      pull(prcp) %>%
      hist(breaks = n, plot = FALSE) -> hist
    hist$counts %>% sum() -> sum
    x <- hist$mids 
    y <- hist$counts/sum 
    data.frame(x, y) -> dens 
    mod <- drm(data = dens, y ~ log(x), fct = EXD.2())
    d <- mod$parmMat[1]
    e <- mod$parmMat[2]
    k <- 1/e/d
    p <- 1 - pexp(log(600)*k, d)
    list <- c(list, p)
    listB <- c(listB, n)
    Year <- c(Year, i)
    dlist <- c(dlist, d)
    elist <- c(elist, e)
  }}
mtr <- data.frame(year = Year,
                    Breaks = listB, 
                      Probability = list, 
                  parameter_d = dlist,
                  parameter_e = elist)

mtr %>%
  mutate(Breaks = as.factor(Breaks)) %>%
  ggplot(aes(x = year, y = Probability, color = Breaks)) + 
    geom_point() + 
    geom_line() + 
    scale_color_brewer(palette = "Blues") +
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        )

g <- df %>%
  ggplot(aes(x = date, y = prcp)) + 
  geom_line(color = "blue") + 
  theme(panel.border = element_blank(),
          panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_rect(fill = "#fafafa"),
          plot.background = element_rect(fill = "white")
        ) + 
  geom_line(aes(x = date, y = rep(1500, length(date))), color = "lightblue") + 
  ylab("Precipitation (in tenth of mm")
ggplotly(g)

data %>%
  filter(prcp >= 1000)

mtr %>%
  filter(Breaks == 21) %>%
  left_join(meth, by = c("year" = "Year")) %>%
  select(year, Probability, Length) %>%
  mutate(days = Probability*Length)